const double mpi = (139.5706/931.494061-0.04)*931.494061;
const double mk = (493.677/931.494061-0.04)*931.494061;
const double mp = (1.00728-0.04)*931.494061;
const double md = (2.01355-0.04)*931.494061;
const double mt = (3.0155-0.04)*931.494061;
const double pi = 3.14159265358979;

void prim(TAxis *xx, const char *title);
inline double intg(double *p){ return fabs(p[0]*sqrt(2.*pi)*p[2])/10.; }

// #define TRK3D

#ifdef TRK3D
string rootfile = "pion3d.root";
#else
string rootfile = "pion.root";
#endif

void fitK(){
  TCanvas *c = new TCanvas("c", "c", 800, 600);
  gStyle->SetOptStat(0);
  TFile *f = new TFile(rootfile.c_str());

  TTree *etf = (TTree *)f->Get("etf");
  etf->Draw("pidm.m()>>hm(1000, -4000., 6000.)", "", "");
  TH1F *hm = (TH1F *)f->Get("hm");
  const int n = hm->GetEntriesFast();
  cout << "hm->GetEntriesFast(): " << n << endl;
  hm->SetTitle("");
  TAxis *xx = hm->GetXaxis();
  xx->SetTitle("m/Z [MeV]"); xx->CenterTitle();
  xx->SetTitleSize(0.05); xx->SetTitleOffset(0.88);
  c->Print("m.pdf");

  // commence the fit //
  TF1 *gpi = new TF1("gmpi", "gaus", 0., 6000.);
  // TF1 *gk = new TF1("gmk", "gaus", 0., 6000.);
  TF1 *gp = new TF1("gmp", "gaus", 0., 6000.);
  TF1 *gd = new TF1("gmd", "gaus", 0., 6000.);
  TF1 *gt = new TF1("gmt", "gaus", 0., 6000.);
  TF1 *bg = new TF1("gbg", "pol1", 0., 6000.);

  gpi->SetParameter(1, mpi); gpi->SetParLimits(1, mpi-100., mpi+100.);
  // gk->SetParameter(1, mk); gk->SetParLimits(1, mk-100., mk+100.);
  gp->SetParameter(1, mp); gp->SetParLimits(1, mp-100., mp+100.);
  gd->SetParameter(1, md); gd->SetParLimits(1, md-100., md+100.);
  gt->SetParameter(1, mt); gt->SetParLimits(1, mt-100., mt+100.);
  bg->SetParameter(0, 0.); bg->SetParameter(1, 0.);

  TF1 *tot = new TF1("tot", "gaus(0)+gaus(3)+gaus(6)+gaus(9)+pol1(12)", 0., 6000.);
  // TF1 *tot = new TF1("tot", "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)+pol1(15)", 0., 6000.);
  tot->SetNpx(6000); tot->SetLineColor(kRed); tot->SetLineWidth(2);
  double par[18]{};
  hm->Fit(gpi, "NR", "", 0., 200.);
  // hm->Fit(gk, "NR+", "", 250., 650.);
  hm->Fit(gp, "NR+", "", 650., 1350.);
  hm->Fit(gd, "NR+", "", 1500., 2500.);
  hm->Fit(gt, "NR+", "", 2200., 3200.);
  hm->Fit(bg, "NR+", "", 0., 6000.);
  gpi->GetParameters(&par[0]);
  // gk->GetParameters(&par[3]);
  gp->GetParameters(&par[3]);
  gd->GetParameters(&par[6]);
  gt->GetParameters(&par[9]);
  bg->GetParameters(&par[12]);
  tot->SetParameters(par);
  tot->SetParLimits(1, 0., 190.);
  tot->SetParLimits(12, 0., 50.);
  // tot->SetParLimits(4, 400., 600.);
  // tot->SetParLimits(12, 0., 50.);
  tot->SetParLimits(0, 0., n/90.);
  // tot->SetParLimits(3, 0., n/180.);
  tot->SetParLimits(3, 0., n/50.);
  tot->SetParLimits(6, 0., n/180.);
  tot->SetParLimits(9, 0., n/180.);
  tot->SetParLimits(1, mpi-100., mpi+100.);
  // tot->SetParLimits(4, mk-200., mk+200.);
  tot->SetParLimits(4, mp-300., mp+300.);
  tot->SetParLimits(7, md-400., md+400.);
  tot->SetParLimits(10, mt-500., mt+500.);
  hm->Fit(tot, "R+", "", 0., 6000.);
  tot->GetParameters(par);
  // for(auto p : par) cout << " " << p;
  // cout << endl;
  gpi->SetParameters(&par[0]); gpi->Draw("same");
  // gk->SetParameters(&par[3]); gk->Draw("same");
  gp->SetParameters(&par[3]); gp->Draw("same");
  gd->SetParameters(&par[6]); gd->Draw("same");
  gt->SetParameters(&par[9]); gt->Draw("same");
  bg->SetParameters(&par[12]); bg->Draw("same");
  gpi->SetLineColor(3); gpi->SetLineStyle(2);
  // gk->SetLineColor(4); gk->SetLineStyle(2);
  gp->SetLineColor(9); gp->SetLineStyle(2);
  gd->SetLineColor(6); gd->SetLineStyle(2);
  gt->SetLineColor(7); gt->SetLineStyle(2);
  bg->SetLineColor(8); bg->SetLineStyle(2);

  // fit the minus-charged particles //
  TF1 *gpim = new TF1("gmpim", "gaus", -4000., 0.);
  TF1 *gkm = new TF1("gmkm", "gaus", -4000., 0.);
  // TF1 *gpm = new TF1("gmpm", "gaus", -4000., 0.);
  // TF1 *gdm = new TF1("gmdm", "gaus", -4000., 0.);
  // TF1 *gtm = new TF1("gmtm", "gaus", -4000., 0.);
  TF1 *bgm = new TF1("gbgm", "pol1", -4000., 0.);

  gpim->SetParameter(1, -mpi); gpim->SetParLimits(1, -mpi-100., -mpi+100.);
  gkm->SetParameter(1, -mk); gkm->SetParLimits(1, -mk-100., -mk+100.);
  // gpm->SetParameter(1, -mp); gpm->SetParLimits(1, -mp-100., -mt+100.);
  bgm->SetParameter(0, 0.); bgm->SetParameter(1, 0.);

  TF1 *totm = new TF1("totm", "gaus(0)+gaus(3)+pol1(6)", -4000., 0.);
  // TF1 *totm = new TF1("totm", "gaus(0)+gaus(3)+gaus(6)+gaus(9)+gaus(12)+pol1(15)", -4000., 0.);
  totm->SetNpx(6000); totm->SetLineColor(kRed); totm->SetLineWidth(2);
  double parm[17]{};
  hm->Fit(gpim, "NR", "", -200., 0.);
  hm->Fit(gkm, "NR+", "", -650., -250.);
  // hm->Fit(gpm, "NR+", "", -1350., -650.);
  // hm->Fit(gdm, "NR+", "", -2500., -1500.);
  // hm->Fit(gtm, "NR+", "", -3200., -2200.);
  hm->Fit(bgm, "NR+", "", -4000., 0.);
  gpim->GetParameters(&parm[0]);
  gkm->GetParameters(&parm[3]);
  // gpm->GetParameters(&parm[6]);
  // gdm->GetParameters(&parm[9]);
  // gtm->GetParameters(&parm[12]);
  bgm->GetParameters(&parm[6]);
  totm->SetParameters(parm);
  totm->SetParLimits(1, -150., 0.);
  totm->SetParLimits(6, 0., 50.);
  totm->SetParLimits(0, 0., n/250.);
  totm->SetParLimits(3, 0., n/380.);
  // totm->SetParLimits(6, 0., n/770.);
  // totm->SetParLimits(9, 0., n/770.);
  // totm->SetParLimits(12, 0., n/770.);
  totm->SetParLimits(1, -mpi-100., -mpi+100.);
  totm->SetParLimits(4, -mk-200., -mk+200.);
  // totm->SetParLimits(7, -mp-300., -mp+300.);
  // totm->SetParLimits(10, -md-400., -md+400.);
  // totm->SetParLimits(13, -mt-500., -mt+500.);
  hm->Fit(totm, "R+", "", -4000., 0.);
  totm->GetParameters(parm);
  gpim->SetParameters(&parm[0]); gpim->Draw("same");
  gkm->SetParameters(&parm[3]); gkm->Draw("same");
  // gpm->SetParameters(&parm[6]); gpm->Draw("same");
  // gdm->SetParameters(&parm[9]); gdm->Draw("same");
  // gtm->SetParameters(&parm[12]); gtm->Draw("same");
  bgm->SetParameters(&parm[6]); bgm->Draw("same");
  gpim->SetLineColor(2); gpim->SetLineStyle(2);
  gkm->SetLineColor(3); gkm->SetLineStyle(2);
  // gpm->SetLineColor(4); gpm->SetLineStyle(2);
  // gdm->SetLineColor(5); gdm->SetLineStyle(2);
  // gtm->SetLineColor(6); gtm->SetLineStyle(2);
  bgm->SetLineColor(7); bgm->SetLineStyle(2);
  cout << "Integrals:" << endl;
  // cout << gpim->Integral(-5000., 5000.) << endl;
  cout << "pi: " << intg(&par[0]) << endl;
  cout << "kaon: " << intg(&par[3]) << endl;
  cout << "p: " << intg(&par[6]) << endl;
  cout << "d: " << intg(&par[9]) << endl;
  cout << "t: " << intg(&par[12]) << endl;
  cout << "-----------------" << endl;
  cout << "pim: " << intg(&parm[0]) << endl;
  cout << "kaonm: " << intg(&parm[3]) << endl;
  // cout << "pm: " << intg(&parm[6]) << endl;
  // cout << "dm: " << intg(&parm[9]) << endl;
  // cout << "tm: " << intg(&parm[12]) << endl;
  cout << "# entries: " << n << " events" << endl;
  int sum = intg(&par[0]) + intg(&par[3]) + intg(&par[6]) + intg(&par[9]) + intg(&par[12]);
  int summ = intg(&parm[0]) + intg(&parm[3]);
  // int summ = intg(&parm[0]) + intg(&parm[3]) + intg(&parm[6]) + intg(&parm[9]) + intg(&parm[12]);
  cout << "tot: " << sum+summ << endl;

  // return;

  TLegend *lg = new TLegend(0.6, 0.4, 0.89, 0.89);
  lg->SetNColumns(2);
  lg->SetBorderSize(0);
  lg->AddEntry(tot, "total", "l");
  lg->AddEntry(gpi, "\\pi^{+}", "l");
  lg->AddEntry(gpim, "\\pi^{-}", "l");
  // lg->AddEntry(gk, "K^{+}", "l");
  lg->AddEntry(gkm, "K^{-}", "l");
  lg->AddEntry(gp, "p", "l");
  // lg->AddEntry(gpm, "\\bar{p}", "l");
  lg->AddEntry(gd, "d", "l");
  // lg->AddEntry(gdm, "\\bar{d}", "l");
  lg->AddEntry(gt, "background", "l");
  // lg->AddEntry(gtm, "\\bar{t}", "l");
  lg->AddEntry(bg, "background", "l");
  // for the minus particles //
  // lg->AddEntry(totm, "total", "l");
  // lg->AddEntry(bgm, "background", "l");
  lg->Draw();
  c->Print("fit.pdf");

  return;

  // x1 v.s. x2 //
  etf->Draw("pidm.x2:pidm.x1>>h1(500, -1200., 1200., 500, -1200., 1200.)", "", "col");
  TH1F *h1 = (TH1F *)f->Get("h1");
  h1->SetTitle("x2 v.s. x1");
  prim(h1->GetXaxis(), "x1 [mm]");
  prim(h1->GetYaxis(), "x2 [mm]");
  c->Print("x12.pdf");

  // pid //
  TF1 *fpionplus = new TF1("fpionplus", "1.+pow(139.5706/x, 2)", 0., 5000.);
  fpionplus->SetNpx(10000);
  TF1 *fpionmius = new TF1("fpionmius", "1.+pow(139.5706/(-x), 2)", -5000., 0.);
  fpionmius->SetNpx(10000);
  TF1 *fkaonplus = new TF1("fkaonplus", "1.+pow(493.677/x, 2)", 0., 5000.);
  fkaonplus->SetNpx(10000);
  TF1 *fproton = new TF1("fproton", "1.+pow(1.00728*931.494061/x, 2)", 0., 5000.);
  fproton->SetNpx(10000);
  TF1 *fdeutron = new TF1("fdeutron", "1.+pow(2.01355*931.494061/x, 2)", 0., 5000.);
  fdeutron->SetNpx(10000);
  TF1 *ftriton = new TF1("ftriton", "1.+pow(3.0155*931.494061/x, 2)", 0., 5000.);
  ftriton->SetNpx(10000);
  // energy spectrum of pion
  // 1. get a cut of pionplus: pionplus
  etf->Draw("pidm.bm2_(0):pidm.poz>>h2(500, -1000., 2000., 500, 0., 20)", "", "col");
  TH2F *h2 = (TH2F *)f->Get("h2");
  h2->SetTitle("\\beta^{-2} v.s. p/Z");
  prim(h2->GetXaxis(), "p/Z [MeV/c]");
  prim(h2->GetYaxis(), "\\beta^{-2}");
  c->Print("pidm0.pdf");
  fpionplus->Draw("same");
  fpionmius->Draw("same");
  fkaonplus->Draw("same");
  fproton->Draw("same");
  fdeutron->Draw("same");
  ftriton->Draw("same");
  c->Print("pidm.pdf");

  c->SetLogy();
  const char *pionCut = "pidm.aoz>0.&&pidm.aoz<0.25"; // pion
  etf->Draw("pidm.pozt>>h3(500, -50., 500.)", pionCut);
  TH1F *h3 = (TH1F *)f->Get("h3");
  h3->SetTitle("poz_{\\perp} of  \\pi^{+}");
  prim(h3->GetXaxis(), "poz_{\\perp} [MeV/c]");
  c->Print("pt_pion.pdf");

  etf->Draw("pidm.pozp>>h3p(500, -50., 500.)", pionCut);
  TH1F *h3p = (TH1F *)f->Get("h3p");
  h3p->SetTitle("poz_{\\parallel} of  \\pi^{+}");
  prim(h3p->GetXaxis(), "poz_{\\parallel} [MeV/c]");
  c->Print("pp_pion.pdf");

  etf->Draw("pidm.beta>>h3b(500, 0., 1.)", pionCut);
  TH1F *h3b = (TH1F *)f->Get("h3b");
  h3b->SetTitle("\\beta of \\pi^{+}");
  prim(h3b->GetXaxis(), "\\beta");
  c->Print("b_pion.pdf");


  const char *pionmCut = "pidm.aoz<0.&&pidm.aoz>-0.25"; // pionm
  etf->Draw("pidm.pozt>>h3m(500, -500., 50.)", pionmCut);
  TH1F *h3m = (TH1F *)f->Get("h3m");
  h3m->SetTitle("poz_{\\perp} of  \\pi^{-}");
  prim(h3m->GetXaxis(), "poz_{\\perp} [MeV/c]");
  c->Print("pt_pionm.pdf");

  etf->Draw("pidm.pozp>>h3mp(500, -500., 50.)", pionmCut);
  TH1F *h3mp = (TH1F *)f->Get("h3mp");
  h3mp->SetTitle("poz_{\\parallel} of  \\pi^{-}");
  prim(h3mp->GetXaxis(), "poz_{\\parallel} [MeV/c]");
  c->Print("pp_pionm.pdf");

  etf->Draw("pidm.beta>>h3mb(500, 0., 1.)", pionmCut);
  TH1F *h3mb = (TH1F *)f->Get("h3mb");
  h3mb->SetTitle("\\beta of \\pi^{-}");
  prim(h3mb->GetXaxis(), "\\beta");
  c->Print("b_pionm.pdf");


  const char *protonCut = "pidm.aoz>0.65&&pidm.aoz<1.35"; // proton
  etf->Draw("pidm.pozt>>h4(500, -50., 500.)", protonCut);
  TH1F *h4 = (TH1F *)f->Get("h4");
  h4->SetTitle("poz_{\\perp} of proton");
  prim(h4->GetXaxis(), "poz_{\\perp} [MeV/c]");
  c->Print("pt_proton.pdf");

  etf->Draw("pidm.pozp>>h4p(500, -50., 2000.)", protonCut);
  TH1F *h4p = (TH1F *)f->Get("h4p");
  h4p->SetTitle("poz_{\\parallel} of proton");
  prim(h4p->GetXaxis(), "poz_{\\parallel} [MeV/c]");
  c->Print("pp_proton.pdf");

  etf->Draw("pidm.beta>>h4b(500, 0., 1.)", protonCut);
  TH1F *h4b = (TH1F *)f->Get("h4b");
  h4b->SetTitle("\\beta of proton");
  prim(h4b->GetXaxis(), "\\beta");
  c->Print("b_proton.pdf");

  c->SetLogy(0);
  etf->Draw("postMagTrackerX.GetFiredStripId()>>h5(33, -1.5, 31.5)");
  TH1F *h5 = (TH1F *)f->Get("h5");
  h5->SetTitle("Fired Strip Id Distribution of TOFWall");
  prim(h5->GetXaxis(), "TOFWall Strip Id");
  c->Print("tofwDist.pdf");

  etf->Draw("mpostMagTrackerX.GetFiredStripId()>>h6(33, -1.5, 31.5)");
  TH1F *h6 = (TH1F *)f->Get("h6");
  h6->SetTitle("Fired Strip Id Distribution of TOFWallm");
  prim(h6->GetXaxis(), "TOFWallm Strip Id");
  c->Print("tofwmDist.pdf");

#ifdef TRK3D
  etf->Draw("postMagTracker3D.GetHitPosY(\"tofw\"):postMagTracker3D.GetHitPosX(\"tofw\")>>"
    "h7(500, -500., 500., 500, -500., 500.)", "", "col");
  TH2F *h7 = (TH2F *)f->Get("h7");
  h7->SetTitle("Beam Position on TOFWall");
  prim(h7->GetXaxis(), "TOFWall X [mm]");
  prim(h7->GetYaxis(), "TOFWall Y [mm]");
  c->Print("tofpos.pdf");

  etf->Draw("pidm.y2():pidm.x2>>h8(500, -800., 800., 500, -400., 400.)", "", "colbox");
  TH2F *h8 = (TH2F *)f->Get("h8");
  h8->SetTitle("Beam Position on TOFWall");
  prim(h8->GetXaxis(), "TOFWall X [mm]");
  prim(h8->GetYaxis(), "TOFWall Y [mm]");
  c->Print("y2x2.pdf");
#endif
} // end function fit



void prim(TAxis *xx, const char *title){
  xx->SetTitle(title); xx->CenterTitle();
  xx->SetTitleSize(0.05); xx->SetTitleOffset(0.8);
} // end function prim